#include "spiral.h"
#include "mask.h"
#include "lines.h"
#include <math.h>

int ArchiSpiral(double l, double Rmin, double g, double w, double o)
{ //draw 1 in 1 out Archimedean spiral with uniform waveguide width and waveguide gap.
  //Archimedean spiral R0+s/2*pi*theta.
  //l: total spiral length
  //Rmin: mininal bent radius
  //g: gap between two waveguide
  //w: waveguide width
  //o: s-bend offset
  
  return ArchiSpiralN(l, Rmin, g, w, o, 1);
  
  /* //1in1out implementation
  int Nh = 10;
  int Ns = 2*Nh; 
  double pi = 4*atan(1.0);
  
  double s = 2*(g + w); //center to center distance
  double astep = 2*pi/Ns; //angle increment
  double Rstep = (s)/Ns; // radius increment
  double Rmin2 = Rmin + o; // take into account the bending offset
  double dy = Rstep/tan(astep*0.5); // dy offset for s-bend
  double dx = 2*(sqrt((Rmin2)*(Rmin2)-(dy*dy*0.25))+(Rmin2)); // dx offset for s-bend
  double R0 = 0.5*(dx+Rstep); // R0 in Archimedean formula
  double a0 = pi - asin(dy*0.5/(Rmin2)); // s-bend turn angle
  
  double l0 = Rmin*a0*2+(R0+R0+(s)*0.5-Rstep)*astep*Nh*0.5;// initial length including s-bend + a half internal loop.
  
  // fprintf(stderr, "%e, %e, %e, %e, %e\n", astep, Rstep, R0, a0, l0);
  
  if (l0 > l)
  {
    fprintf(stderr, "Error: length should be larger than %e for Rmin=%e\n", l0, Rmin);
    return 0;
  }
  
  double Rt = R0;
  int step = 0;
  
	while (l0 < l)  {
		l0 += Rt*astep + (Rt+0.5*s)*astep;
    Rt += Rstep;
    step++;
	}
  Rt -= Rstep;
  step--;
  
  double adif = astep - (l0-l)/(2*Rt+0.5*s);
  fprintf(stderr, "%e, %e, %d\n", adif, Rt, step);
  
  //main spiral
  push_cp();
  Rt = Rt + 0.25*s;
  cw2(Rt, adif, w, s*0.5-w);
  for(int i=0; i<step; i++)
  {
    Rt -= Rstep;
    cw2(Rt, astep, w, s*0.5-w);
  }
  
  //additional half
  push_cp();
  offset(-0.25*s);
  Rt = Rt + 0.25*s;
  for(int i=0; i<Nh; i++)
  {
    Rt -= Rstep;
    cw(Rt, astep, w);
  }
  
  // s-bend
  pop_cp();
  offset(0.25*s+o);
  
  cw(Rmin, a0, w);
  offset(-2*o);
  cw(Rmin, -a0, w);
  pop_cp();
  rotate(RAD(180));
  return 1;
  */ ////1in1out implementation
  
}

int ArchiSpiralN(double l, double Rmin, double g, double w, double o, int Nline)
{ //draw #Nline in #Nline out Archimedean spiral with uniform waveguide width and waveguide gap.
  //Archimedean spiral R0+s/2*pi*theta.
  //l: total spiral length
  //Rmin: mininal bent radius
  //g: gap between two waveguide
  //w: waveguide width
  //o: s-bend offset
  //Nline: the input waveguide number. (2Nline = input + output)

  
  int Nh = 10;// segments of the pi degree curve.
  int Ns = 2*Nh; //segments of the 2*pi degree curve.
  double pi = 4*atan(1.0); //get pi
  
  double s = 2*(g + w)*Nline; //center to center distance, s in Archimedean formula.
  
  double astep = 2*pi/Ns; //angle increment
  double Rstep = (s)/Ns; // radius increment
  double Rmin2 = Rmin + o + ((Nline/2) + (Nline%2))*0.5*(g+w) ; // take into account the bending offset, (the middle radius of the line cluster)
  double dy = Rstep/tan(astep*0.5); // dy offset for s-bend for the middle
  double dx = 2*(sqrt((Rmin2)*(Rmin2)-(dy*dy*0.25))+(Rmin2)); // dx offset for s-bend for the middle
  double R0 = 0.5*(dx+Rstep); // R0 in Archimedean formula
  double a0 = pi - asin(dy*0.5/(Rmin2)); // s-bend turn angle
  
  double l0 = Rmin*a0*2+(R0+R0+(s)*0.5-Rstep)*astep*Nh*0.5;// initial length including s-bend + a half internal loop.
  
  // fprintf(stderr, "%e, %e, %e, %e, %e\n", astep, Rstep, R0, a0, l0);
  
  if (l0 > l)//
  {
    fprintf(stderr, "Error: length should be larger than %e for Rmin=%e\n", l0, Rmin);
    return 0;
  }
  
  double Rt = R0;//temp radius
  int step = 0;//turn steps
  
	while (l0 < l)  {
		l0 += Rt*astep + (Rt+0.5*s)*astep;
    Rt += Rstep;
    step++;
	}
  Rt -= Rstep;
  step--;
  
  double adif = astep - (l0-l)/(2*Rt+0.5*s);// last turn angle (<= astep)
  //fprintf(stderr, "%e, %e, %d\n", adif, Rt, step);
  
  push_cp();
  Rt = Rt + 0.25*s;
  cwN(Rt, adif, w, g, Nline*2);
  for(int i=0; i<step; i++)
  {
    Rt -= Rstep;
    cwN(Rt, astep, w, g, Nline*2);//
  }
  push_cp();
  offset(-0.25*s);
  Rt = Rt + 0.25*s;
  for(int i=0; i<Nh; i++)
  {
    Rt -= Rstep;
    cwN(Rt, astep, w, g, Nline);
  }
  
  pop_cp();
  offset(0.25*s+o);
  
  cwN(Rmin2-o, a0, w, g, Nline);
  offset(-2*o);//s-bend offset
  cwN(-(Rmin2-o), -a0, w, g, Nline);
  
  pop_cp();
  rotate(RAD(180));//final cp stopped at the middle of the input but rotate the direction by 180C
  return 1;
}


